#'###########################################
#'                                          #
#'  ------------------------------------    #
#'  |     COIN     __                  |    #
#'  |            >(' )                 |    #
#'  |              )/                  |    #
#'  |             /  `----/            |    #
#'  |             \  ~=- /             |    #
#'  |   ~^~^~^~^~^~^~^~^~^~^~^~^~^~^   |    #
#'  ------------------------------------    #
#'                                          #
#'###########################################


HRcompare.mobile2 <- function(errors, gps, iterations = 100, plot = FALSE, resultsPath=".", write.results=FALSE, file.name="",
		mcp.percent=100, kernel.percent=99, kernel.h="href", kernel.cellsize=250, kernel.extent=0.5) {
	print("Comparing Home ranges")
	
	errors <- errors[, list(argLong, argLat, LC, time, deployment)]
	
	gps[, time := as.POSIXct(time, format="%d/%m/%Y %H:%M:%S", tz="EST5EDT")]
	setkey(errors, deployment, LC)
	setkey(gps, deployment, time)
	
	
	#' Comparaison des domaines vitaux par deploiement
	res <- errors[, compareHomeRanges2(.SD, gps, iterations, .BY[[1]], mcp.percent, kernel.percent, kernel.h, kernel.cellsize, kernel.extent), by=deployment]
	
	if (write.results){
		filename <- paste("HR_compare_", file.name, ".csv", sep="")
		write.csv(res, file.path(resultsPath, filename))
	}
	return (res)
}

compareHomeRanges2 <- function (argoslocs, gps, iterations, depl, mcp.percent, kernel.percent, kernel.h, 
		kernel.cellsize, kernel.extent) {
	print(paste("Computing for deployment", depl))
	
	n <- nrow(argoslocs)
	ncols <- length(HRCOMPARE_COLNAME) / 2
	
	#' si on a moins de 5 localisations, on ne peut pas faire de DV
	if (n < 5) {
		res <- as.list(as.double(c(n, rep(NA, length(HRCOMPARE_COLNAME) - 1))))
		names(res) <- HRCOMPARE_COLNAME
		return (res)
	}
	
	#' on ne garde que les donnees gps dans le creneau horaire argos
	gpsdata <- gps[deployment == depl]
	gpsdata <- gpsdata[time < max(argoslocs$time) & time > min(argoslocs$time)]
	
	#' Calcul du domaine vital argos via MCP. Inutile de le faire plus d'une fois
	argos_sp <- SpatialPoints(argoslocs[, list(argLong, argLat)])
	argos.mcp <- mcp(argos_sp, percent=mcp.percent)
	argosArea.mcp <- argos.mcp@polygons[[1]]@area
	
	#' Calcul du domaine vital argos via kernel
	kernel.grid <- makegrid.kernel(argoslocs, gpsdata, kernel.cellsize, kernel.extent)
	argos.kernel <- kernelUD(argos_sp, h=kernel.h, grid=kernel.grid)
	
	
	#' Multithreading si possible pour faire les iterations
	#if (exists("ARGOS_CLUSTER")) {
	# rep <- iterations / ARGOS_CLUSTER_CORES
	#	tmp <- clusterCall(ARGOS_CLUSTER, compare.gps, gpsdata, rep, ncols, n, mcp.percent, kernel.percent,
	#			kernel.h, kernel.grid, argos.mcp, argos.kernel, argosArea.mcp)
	#	tmp <- do.call(rbind, tmp)
	#} else {
	tmp <- compare.gps2(gpsdata, iterations, ncols, n, mcp.percent, kernel.percent, kernel.h, kernel.grid, argos.mcp, argos.kernel, argosArea.mcp, argos_sp@coords)
	#}
	#' formattage des resultats
	
	res <- c(round(apply(tmp, 2, mean, na.rm=TRUE), 2), round(apply(tmp, 2, sd, na.rm=TRUE), 2))
	names(res) <- HRCOMPARE_COLNAME
	res <- as.list(res)
	return (res)
}



compare.gps2 <- function(gpsdata, rep, res.ncol, n, mcp.percent, kernel.percent, kernel.h,
		kernel.grid, argos.mcp, argos.kernel, argosArea.mcp, argos.coords) {
	tmp <- matrix(0, rep, res.ncol)
	for (j in 1:rep) {
		print(paste("Iteration no: ", j))
		
		#' Sous echantillonnage d'un nombre de points gps egal au nombre de localisations argos
		gpsindexes <- sample(nrow(gpsdata), n)
		gpslocs <- gpsdata[gpsindexes, list(longitude, latitude)]
		gps_sp <- SpatialPoints(gpslocs)
		
		
		#'#
		#'# MCP 
		#'# 
		if (mcp.percent > 0) {
			gps.mcp <- mcp(gps_sp, percent=mcp.percent)
			
			#' Aire
			gpsArea.mcp <- gps.mcp@polygons[[1]]@area
			
			#' ratio de taille argos/gps
			sizeRatio.mcp <- argosArea.mcp / gpsArea.mcp
			
			#' Surface d'overlap nulle si pas d'intersection
			intersection.mcp <- gIntersection(gps.mcp, argos.mcp)
			if (is.null(intersection.mcp)){
				overlapArea.mcp <- 0
			} else {
				overlapArea.mcp <- intersection.mcp@polygons[[1]]@area
			}
			
			#' Index de chevauchement de Dahl (2005)
			overlapIndex.mcp <- (overlapArea.mcp / (argosArea.mcp + gpsArea.mcp - overlapArea.mcp)) * 100
			
			#' Pourcentage de points dans le polygone
			ispip.mcp <- point.in.polygon(argos.coords[, 1], argos.coords[, 2], 
					gps.mcp@polygons[[1]]@Polygons[[1]]@coords[, 1], 
					gps.mcp@polygons[[1]]@Polygons[[1]]@coords[, 2])
			prop.in.gps.mcp <- sum(ispip.mcp) / n * 100
		}
		
		#'#
		#'# Kernel
		#'# 
		if (kernel.percent > 0) {
			#' DV avec Kernel pour GPS
			gps.kernel <- kernelUD(gps_sp, h=kernel.h, grid=kernel.grid)
			
			#' Creation d'une fausse classe estUDm pour calculer l'overlap
			kernelHR <- list(argos=argos.kernel, gps=gps.kernel)
			class(kernelHR) <- "estUDm"
			
			#' Aires kernel
			argos.kernel.vertices <- getverticeshr(argos.kernel, percent=kernel.percent, unout="m2")
			gps.kernel.vertices <- getverticeshr(gps.kernel, percent=kernel.percent, unout="m2")
			print(as.data.frame(argos.kernel.vertices))
			
			argosArea.kernel <- argos.kernel.vertices@data$area
			gpsArea.kernel <- gps.kernel.vertices@data$area
			
			#' ratio de taille argos/gps
			sizeRatio.kernel <- argosArea.kernel / gpsArea.kernel
			
			#' Surface d'overlap nulle si pas d'intersection
			intersection.kernel  <- gIntersection(gps.kernel.vertices, argos.kernel.vertices)
			if (is.null(intersection.kernel)){
				overlapArea.kernel <- 0
			} else {
				overlapArea.kernel <- intersection.kernel@polygons[[1]]@area
			}
			
			#' Index de chevauchement de Dahl (2005)
			overlapIndex.kernel <- (overlapArea.kernel / (argosArea.kernel + gpsArea.kernel - overlapArea.kernel)) * 100
			
			#' Pourcentage de points dans le polygone
			ispip.kernel <- point.in.polygon(argos.coords[, 1], argos.coords[, 2], 
					gps.kernel.vertices@polygons[[1]]@Polygons[[1]]@coords[, 1], 
					gps.kernel.vertices@polygons[[1]]@Polygons[[1]]@coords[, 2])
			prop.in.gps.kernel <- sum(ispip.kernel) / n * 100
		}
		## image(argos.kernel)
		## plot(argos.kernel.vertices, border="brown", add=T)
		## points(argos.coords[, 1], argos.coords[, 2])
		## points(gpslocs, col="green")
		## plot(gps.mcp, add=T)
		## plot(gps.kernel.vertices, border="blue", add=T)
		## plot(getverticeshr(gps.kernel, percent=50), border="red", add=T)
		## plot(getverticeshr(argos.kernel, percent=50), border="pink", add=T)
		
		
		#str(kernel.poly)
		#print(ispip)
		
		tmp[j,] <- c(n, sizeRatio.mcp, overlapIndex.mcp, sizeRatio.kernel, overlapIndex.kernel, argosArea.mcp, gpsArea.mcp, overlapArea.mcp, 
				argosArea.kernel, gpsArea.kernel, prop.in.gps.mcp, prop.in.gps.kernel)
	}
	return (tmp)
}

